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' We report precise simulations of (j) theory in the Ising hmit. A recent technique 



to stochasticahy evaluate the all-order strong coupling expansion is combined with 
exact identities in the closely related Aizenman random current representation. In 
D I this way estimates of the renormalized coupling close to the continuum limit become 

possible with unprecedented precision and yet low CPU cost. As a sample application 
. we present results for the unbroken phase of the Ising model in dimensions 3, 4 and 5 

o : 

, and investigate the question of triviality by studying a finite size scaling continuum 

' limit. 

(N 
O , 

^ . The 0^ theory of a real scalar field is the starting point of many textbooks on quantum 

field theory. It also plays a phenomenological role as an extremely simplified model of the 



^ ■ Higgs sector of the standard model. In Q], building on 2|], it was discussed that the theory 

H ; 

. , is trivial in dimension four, meaning that no interaction can exist in the continuum limit. 

For applications this means, that the effective theory with an unremovable cutoff in place 
has only a limited energy domain of validity. This important property is rigorously known 
to hold above four dimensions j^, Q], but in D = 4 we have to rely so-far on numerical 
checks. In this letter we report on the discovery of a new numerical strategy and algorithm, 
which enhances the precision of such checks by orders of magnitude as the continuum limit is 
approached. Also novel is the use of a finite volume renormalization scheme in this context. 
Due to both improvements we can progress much deeper into the universal scaling region 
than in previous computations like We find a 'borderline' agreement with standard 
perturbation theory for the cutoff dependence of the interaction strength which calls for 



2 



more study. 

We here consider the Ising hmit of (j)^ on a D- dimensional hypercubic lattice of extent 
L and lattice spacing a in all directions. We mostly use lattice units a = 1 from here on, 
but occasionally re-introduce a to emphasize cutoff dependencies. Physical information is 
extracted via n-point correlation functions 

with 

Zixi,X2, . . . , x„) = 2-^ 5^ <^>(yhixi)six2) ■■■s (2) 

and the volume V = L^. We sum over all Ising configurations s{x) = ±1 and Z{^) is the 
proper partition function with no field insertions. On our finite lattice Z{.) is analytic for 
all values of /3. We parametrize the strong coupling expansion in f3 by summing in addition 
to s over an integer link field k{l) = 0, 1, . . . , cxd 

X2, . . . , Xn) = 2-^ Yl ^[^] n [six)s{y)]''^'\s{xi)six2) ■■■s (3) 

s,k l={xy) 

with the multiple Poisson weight 

-W = nf(^. (4) 

For each k the s sum can now be performed and leaves behind a constraint for fc, 

X2, . . . , Xn) = Yw[k]5Q[k],X- (5) 
k 

The Kronecker 6 enforces the coincidence of two sets. The source set Q[k] consists of the 
sites surrounded by an odd total number of k{l), 



g[A;] = <jx| ^ = l(mod2) I. (6) 

I l,dl5x ) 

The insertion set X coincides with (xi, X2, . . . , x„) if they are mutually different, but is more 
generally given by 

X=L|^<5,,,, = l(mod2)l. (7) 

n n " 

Michael Aizenman [3| , [6| has used the above representation of the Ising model to obtain 
rigorous correlation inequalities. He calls {k{l)} random currents and the sets Q = X defects 
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or external sources. Among his results the following is of interest here as it can be turned 
into an efficient numerical algorithm. From proposition 5.1 in the identity 

Zc{xi,X2,Xs,X4) = -2'^w[k]w[k']5Qik],x^2^Q[k'],Xi4X{xi,X3;k + k') (8) 

k,k' 

with the connected part of of Z 

Zc{Xi,X2, Xs, Xi) = Z{xi,X2, X3, X4) - Z {xi, X2) Z {X3, X4) 

- Z{xi,X3)Z{x2,Xi) - Z{xi,Xi)Z{x2,X3) (9) 

follows. The sets X12 and X34 are formed as in ([7]) but from only two points each. The 
cluster incidence function A" G {0, 1} is one if Xi, X3 are in the same bond percolation cluster 
built by bonds that are active on links where k{l) + k'{l) > holds in the doubled random 
current system. Note that the pairs {xi,X2} and {x3,X4,} are connected automatically for 
k, k' that contribute. 

In the Monte Carlo community is has recently been found 9] that it is both 

possible and advantageous to simulate the untruncated strong coupling expansion instead 
of the original path integral (or sum) over fields. In a simple variant of the worm algorithm 

one simulates the ensemble corresponding to the partition function 

Z=J2Mk]5Q[kix.. (10) 

u,v,k 

with a corresponding definition of expectation values of observables {{0[u,v, k])) . In the 
sum u, V run over all lattice sites. 

To simulate this ensemble our elementary update step is as follows: 

• pick at random one of the 2D links emanating from u and call it / = (uu) 

• assign a new value k{l) to this link with probability = exp[— /3) (3^ / k\ 

• if k — k is odd, move u ^ u, otherwise leave u unchanged 

If we alternate these steps with similar ones for v we have a correct algorithm for ( |TOl) . 
Ergodicity may be shown by steps deforming an arbitrary configuration to the trivial one. 
This local heatbath has proved to be slightly superior to Metropolis proposals with k{l) — >■ 
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k{l) ± 1. It is not difficult to show that the Ising two-point function is now given by the 
ratio of histograms 

(^(^)^(0)) = %^ (11) 

where we have used translation invariance. This implies in particular that the susceptibility 
is given by 

X2 = E(^(^)^(0)) = [(«'^))]"'- (12) 

X 

To now make use of (El) for the connected four point susceptibility 

X4 = Y1 {six)siy)siz)sm^ (13) 

x,y,z 

with subtractions as in ([9]) all we have to do is simulate two independent replica of (|TOl) and 
sum over all Xj in ([8]) to arrive at 

_ ^ ^ J{Xiu,u',k + k'))) ^ ^^^^ 

{{Su,vSu' ,v')) 



In total we thus have derived 

1 X4 



2{{X{u,u\k + k'))) . (15) 



The right hand side is obviously bounded between and 2. In particular the lower bound 
corresponds to the Lebowitz inequality. Our estimator reflects this property manifestly, the 
subtraction of disconnected parts has been achieved analytically. We expect this to lead 
to a superior precision for compared to conventional Monte Carlo procedures since they 
involve substantial numerical cancellations here with the correspondingly enhanced relative 
errors. 

The previous expression is strongly reminiscent of a standard definition of a dimensionless 
universal renormalized coupling constant in 0^ theory including the Ising limit with its 
infinite bare coupling. It is given by 

9r = —, — T^"^ (16) 

where m is a the renormalized mass. 

Often, for example in , the mass is defined in terms of the two point function in an 

infinite volume at vanishing momentum. We substitute this by a definition using the two 



5 



smallest possible momenta in a periodic volume as in [8']. The two point function (fTT!) in 
momentum space may be measured by 

~ {{cos{p{u - v)))) 

G{p) = TT^^T • (17) 

Our definition of a renormalized mass m is 



m^+pi G{0) 
where we use the smallest momentum 



{{cosip^iu-v)))) (18) 



= (27r/L,0,0, ...,0), ap* = 2sin(7ra/L) (19) 



and average over its D possible directions. The rationale of the definition ([T6|) is that it 
is a dimensionless ratio with the same number of fields in the numerator and denominator 
and hence it is expected to have a universal continuum limit. As it vanishes for Gaussian 
theories it is a measure of the interaction strength. We are thus led to the definition 

g = 2{{X{u,u',k + k'))) X z", z = mL. (20) 

Combining triviality with finite size scaling we investigate the proposition that the con- 
tinuum limit at fixed z forces \ 0. As for other questions on non-perturbative ultraviolet 
renormalization we find it advantageous to employ a finite volume renormal- 

ization scheme also here. We shall perform a sequence of simulations of growing L = L/a 
where we tune (3 such as to maintain a fixed value z = 2. The advantage of keeping L finite 
and not too large in physical length units is that for the manageable values of L/a we 
expect to be closer to the universal continuum limit. If the theory is trivial we should find 
(jf ^ as L/a ^ oo. This is expected j^, jst, jj] for D > 4, and likely, although only at a 
logarithmic rate, for D = A. 

To get confidence in our algorithmic implementation we first reproduced the results of 
cluster simulations in |5| within errors. For the 20"^ lattice our error in g for a comparable 
number of Flops is about 12 times smaller. 

In the table we compile our data. Each line corresponds to 10^ iterations with V link 
updates in each of the two replica. The cost is dominated by the runs L* = 4, L = 32 and 
L) = 5, L = 16 with 240 hours each on a single PC. We refer to a code running under matlab 



but importing random numbers from the ranlux generator 



13| in C Q . We have insisted 
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D 


L/a 


13 


z 




dX/dz 


X{z = 2) 


4 


8 


0.148320 


1.9981(27) 


0.39235(96) 


-0.3200(14) 


0.39175(63) 


4 


10 


0.148748 


1.9949(26) 


0.37256(92) 


-0.3193(14) 


0.37093(62) 


4 


12 


0.148996 


1.9992(26) 


0.35493(91) 


-0.3165(15) 


0.35469(60) 


4 


16 


0.149270 


1.9988(25) 


0.33161(91) 


-0.3129(16) 


0.33125(58) 


4 


22 


0.149449 


2.0085(24) 


0.30831(86) 


-0.3030(16) 


0.31088(57) 


4 


32 


0.149571 


1.9956(24) 


0.29028(83) 


-0.2993(20) 


0.28896(55) 


3 


8 


0.217350 


1.9946(36) 


0.59387(100) 


-0.2929(13) 


0.59228(76) 


3 


10 


0.218560 


1.9942(37) 


0.58634(102) 


-0.2950(13) 


0.58463(76) 


3 


16 


0.220153 


2.0047(37) 


0.57240(109) 


-0.3023(14) 


0.57382(77) 


3 


32 


0.221143 


2.0032(39) 


0.56338(118) 


-0.3076(17) 


0.56435(76) 


5 


8 


0.113052 


2.0041(18) 


0.19126(59) 


-0.2390(13) 


0.19223(45) 


5 


10 


0.113336 


1.9993(16) 


0.16037(53) 


-0.2170(13) 


0.16022(41) 


5 


12 


0.113503 


1.9937(15) 


0.13884(47) 


-0.1957(12) 


0.13760(38) 


5 


16 


0.113674 


1.9918(13) 


0.10944(39) 


-0.1656(12) 


0.10809(33) 



TABLE I: Simulation results for -0 = 3,4, 5. 



on luxury level 2, but this part still accounts for only 4 % of the CPU time. The main code 
could clearly be accelerated substantially. Derivatives of X and z with respect to ln/3 can 
be measured as connected correlations with J2i ^(0 their quotient yields our estimate 
for dX/dz. In the end it emerges as a certain nonlinear function of primarily measured 
observables and it as well as all other errors is estimated by the tools provided in [l^. We 
take its measured value in each data set as a fixed constant and then form, now with X and 
z as functions of primary quantities, the combination 

X{z = 2) = X + {2- z)dX/dz. (21) 

A look at the table shows that this is a tiny but sometimes significant correction. We 
can however safely neglect the error of dX / dz and higher terms of the Taylor expansion. It 
came as a pleasant surprise that even where no systematic correction is needed the statistical 
fluctuations in this combination partially cancel and thus reduce the error. This saves more 
than another factor two in run-time. The compensation is actually plausible: when sampled 
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graphs are 'bigger' than average X goes up but the mass goes down. In the table we note 
that relative errors are practically independent of L/a, which means complete absence of 
critical slowing down. Nevertheless the (short) autocorrelations do have to be taken into 
account when errors are determined. 



12 z = 2 
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L/a = 32 22 16 12 10 8 

Q 1/1 \ I I \ I \ L 

0.05 0.1 

a/L 



FIG. 1: All measured couplings with fits (full and dashed lines) discussed in the text. Errorbars 
are barely visible within the symbols. 

In figure [1] all coupling data are plotted. Fits (full and dashed lines) all have acceptable 
X^- The D = 3 values are almost cutoff independent and approach a finite continuum value. 
The fit is A + B{a/ LY with the corrections to scaling exponent [16] uj = 0.85. Due to the 
flatness, u can vary over a wide range including u = 1. For D = 5 the full line is A + Ba/L 
omitting the L = 8 lattice. The dashed fit is Aa/L + B{a/L)^. A 32^ simulation would be 
of interest to better verify triviality for this case. The D = A data show more curvature and 
the dotted lines just connect the data points. A naive 'eyeball' extrapolation in this plot 
against a/L would probably arrive at a nonzero value while on theoretical grounds we shall 
argue now for the dash-dotted line extrapolating to the origin at a vertical slope. 
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0.25 0.3 0.35 0.4 0.45 0.5 

l/ln(L/a) 



FIG. 2: Numerical couplings together with the renormalization group evolution starting from the 
finest lattice. Dashed lines include lattice artifacts of the /J-function. 

To confront the AD data with the theory we plot against [ln(L/a)] ^ in figure [2l 

We solve the Callan Symanzik equation for the cutoff dependence of the coupling starting 
from L = 32. For this plot we have changed to the coupling 

~g = g[l + {za/Lf/8]-^ = g + 0{a') (22) 

differing by small cutoff effects only. For it in contrast to g there are no tree level artifacts 
in the /3-function and also the cutoff corrections of the one and two loop terms are more 
uniform. We have worked out the lattice perturbation theory for our scheme up to two 
loops and could therefore, by relating to also obtain the three loop term (without cutoff 



effects). Details will be reported elsewhere 



171 ] . We here draw the following conclusions: The 



one loop result is accurate to a few percent for the scale changes considered. For instance, it 
accounts for 97% of the change L/a = 32 — 16. The two loop term has a reasonable relative 
size but the wrong sign. The three loop term is the first one that is scheme dependent and 
hence depends on z. It is very large for z = 2 which suggests that renormalized perturbation 
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theory as an asymptotic expansion here fails to improve the leading order. It rather is at 
its limit with only the one loop approximation being numerically accurate at the percent 
level. Nonetheless it seems convincing to now trust the one loop approximation for a < L/32 
(for z = 2) which implies a vanishing g in the continuum limit. On the way to it, also the 
higher loop terms should eventually cooperate to improve the approximation. The dash- 
dotted curve in figure [1] shows the one loop evolution continued. We plan a more detailed 
discussion of the perturbative series and data for other z values in [l^] . 




1.26 



strong coupling order x10^ 

FIG. 3: The order in (3 of the diagrams sampled for L = 32, D = 4 (n = u only). 

We end on a more technical theme concerning the strong coupling simulation. The order 
of the diagrams that the algorithm has picked to be important for our physics is shown in 
figure [31 The peak is at the order of the correlation volume = {am) ~^ and the width 
seems to be controlled by [Ljaf' . 
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